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We present for the first time time-dependent density-matrix renormalization-group simulations (t- 
DMRG) at finite temperatures. It is demonstrated how a combination of finite-temperature i-DMRG 
and time-series prediction allows for an easy and very accurate calculation of spectral functions in 
one-dimensional quantum systems, irrespective of their statistics, for arbitrary temperatures. This 
is illustrated with spin structure factors of XX and XXX spin-^ chains. For the XX model we can 



compare against an exact solution and for the XXX model (Heisenberg antiferromagnet) 
Bethe Ansatz solution and quantum Monte Carlo data. 
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PACS numbers: 75.10.Pq, 75.40.Mg, 78.47.-p 



I. INTRODUCTION 

Strongly correlated quantum systems continue to 
pose central challenges in theoretical condensed-matter 
physics. In the case of one-dimensional systems, we now 
have a full range of techniques to address static and 
dynamic ground-state properties. However, condensed- 
matter experiments typically work at finite temperatures 
that cannot be simply approximated by the ground-state 
physics, and the low-temperature physics of such systems 
is of high interest of its own. Due to the homogeneity 
of the systems under study in space and time, the ex- 
perimental responses are best represented in momentum- 
frequency space, as, for example, by spin structure fac- 
tors 



■IkC 



dte^(S^t)S^O)}, (1) 



where the average is taken with respect to some finitc- 
T density operator and [S^,S^,] = i8wt^ vl Sj . We will 
refer to functions of this form as spectral functions in 
the following; the proposed method does not depend on 
the specifics. Theoretical approaches to calculate such 
spectral functions with high accuracy are very limited. 
On the other hand, experimental progress makes such 
calculations very timely: For example, due to small neu- 
tron flux, neutron-scattering techniques were in the past 
essentially limited to finding scattering maxima. Now, 
they have advanced to a degree that both the amplitude 
and the lineshape, which contain important information 
on the many-body physics involved, can be investigated 
with high accuracy (e.g. Refs.[l| and[H). 

For finite-temperature quantum Monte Carlo (QMC) 
calculations (e.g. positive-definite path integral 3 -^ or 
stochastic series expansion^ representation), spectral 
functions have to be extracted by analytic continuation 
from imaginary-time results^, which is ill-conditioned 
and numerically challenging. In the context of 



DMRG 7 -^, early approaches for finite-temperature spec- 
tral functions have built on transfer-matrix DMR C 10 ' 11 , 
which gives thermodynamics for homogeneous infinite 
systems at high precision. In a first step, analytic con- 
tinuation techniques as in QMC were employed^ 2 -. In 
a second development, transfer matrices were evolved 
in real-time to access autocorrelation function al 14 ; the 
time scales reachable limit resolution in frequency space. 
While this approach would also be amenable to the 
prediction techniques used in the following to circum- 
vent this issue, the determination of momentum-space 
correlators will be particularly easy with the i-DMRG 
mcthod o 15 ! 16 ' 17 we use. In an unrelated development, 
a ground-state DMRG-like technique has been pro- 
posed that invokes a random-sampling and averaging 
procedure^. Its feasibility for dynamical quantities has 
only been tested for very small systems; it is also limited 
to low temperatures. 

In this paper, we present an application of i-DMRG 
at finite temperatures. It is shown how a combination 
of this technique and the linear prediction metho d 19 i 20 
allows one to obtain very accurate results for spectral 
functions with frequency and momentum dependence in 
the entire temperature range. While the demonstration 
focuses on a few particular cases, it will be obvious that 
the approach generalizes straightforwardly. 



II. METHOD 

A. DMRG at finite temperature 

In recent years, DMRG 7 .^ has emerged as a key method 
for the ground-state physics of strongly correlated one- 
dimensional quantum systems^. The time-evolution of 
pure states in one-dimensional strongly correlated spin, 
bosonic or fermionic models can now be simulated by t- 
DMRGJ&i&ii over times 10-100 times longer than the 
typical inverse energy scale of the problem (e.g. inverse 
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FIG. 1: (Color online.) The purification of a density operator 
p on a Hilbert space TL is a pure state \p) on an enlarged 
Hilbert space TL (g A such that p — tr^|p)(p|. The infinite 
temperature state po oc 1 is obtained by preparing each pair 
of a physical site and the corresponding auxiliary site in a 
maximally entangled state. 



very high accuracy, real-time simulations arc plagued by 
the same limitations as at T = 0: the propagation of exci- 
tations through the system leads to entanglement growth 
in the purified state. As entanglement entropy is related 
roughly exponentially to DMRG resources, this strongly 
limits achievable simulation times, or inversely the w- 
resolution for spectral functions, as those are derived by 
Fourier transformation from the real-time data. In or- 
der to circumvent this limitation at very low numerical 
cost, we adapt a linear prediction tech niq u o 19 i 20 already 
successfully employed at T = in Ref . I27M281 . 



hopping energy). It was quickly demonstrated 21 - that 
this method can be easily extended to the simulation of 
the static and dynamic behaviors at T > by using the 
purificatio n 22 i 23 i 24 of the density operator: Any density 
operator p of some physical system TL can be encoded by 
a pure state of a combined physical and ancillary system, 
\p) ETC® A, such that the density matrix is retained by 
tracing out A, 

p = tv A \p){p\. (2) 

The ancillary system A can be taken as a copy of the 
physical state space TL; Fig. [TJ In the case of interest, 
(unnormalized) thermal density operators pp = e~ l3H ', 
the corresponding purification can be constructed by an 
imaginary-time evolution starting from the purification 
of the (trivial) infinite-T ((3 = 0) density operator p = 1. 
A possible purification for this (3 = ensemble is 

\po) = ®f =1 \po,e) with \p 0i t) = ^2 We) ® We)', (3) 

where {|<7^}} and {|o^)'} denote the bases of the physical 
state space of site I and its associated ancillary state 
space, respectively. With this, we have po = Tr A |po) {po I ■ 
Finite temperatures [3 > are reached by imaginary-time 
evolution 

Pp = e~^ = Tr A \pp)(pp\ with \ P p) = e~^/ 2 \p ). 

Proper normalization is restored by imposing (pp\pp) = 
1. The mixed state pp can then be evolved in time as 

\pp{t)) = e- i6t \p (O)) and fls (t) = Tr^ \p (t))(pp(t)\. 

As a product state, the initial (3 = purification \po) 
is uncorrelated, and hence, for the DMRG simulation, 
can be expressed exactly with block Hilbert spaces of di- 
mension m = 1. Imaginary-time evolution will introduce 
correlations, requiring one to increase m. For the evalu- 
ation of expectation values, both physical and ancillary 
degrees of freedom are traced over. As an example, take 

Tr H (S?(t)S?(0)pp) = Tr H ® A (S?(t)S?(0)\pp)(pp\) 
= (pp\e i6t Sle- i6t S?\pp). (4) 

While this approach has been found to yield thermody- 
namic quantities^ and static correlators^ at T > to 



B. Linear prediction 

For a time series of complex data Xo, x\, . . . , xjy at 
equidistant points in time t n = n ■ At (and i b s := NAt) 
one makes a prediction of xjy+i, ■ ■ ■■ In our case, 

t bs is the (physical) time where the t-DMRG calculation 
was stopped - usually because the computational cost to 
simulate further with a given accuracy has become too 
high. For the data points beyond t = t bs! linear predic- 
tion makes the ansatz 

v 

X n ^ ^ Q,jX n — j. (5) 

1=1 

The (predicted) value x n at time step n is assumed to be a 
linear combination of p previous values . . . , x n - p }. 

Once the coefficients a; arc determined from the known 
data, the ansatz is used to calculate (an approximation 
of) all x n with n > N. 

The coefficients cii are determined by minimizing the 
least-square error in the predictions over a subinterval 
(tabs — ifitjiobs] of the known data, i.e. we want to mini- 
mize 

E = X! \x n - x n \ 2 /w n (6) 

where A/fit is the fitting interval A/fit = {n\t n e (£ bs — 
ifiti^obs]}; an d w n is some weighting function; we choose 
w n = 1 for our simulations. We found tfn = t bs/2 to 
be a robust choice. It is a compromise between choosing 
tat small enough to eliminate spurious short-time behav- 
ior from the true long-time behavior and choosing tfi t 
large to have a good statistics for the fit and allow for a 
large number p of coefficients a,i in the ansatz Min- 
imization of the error E, Eq. ©, with respect to the 
coefficients a, yields the system of linear equations 

Ra = -r, (7) 

where R and r are the autocorrelations 

Rji ^ X n _jX n — i I W n , Tj — ^ X n _jX n J W n . 

Equation ([7]) is solved by a = — i? _1 r. For positive w n , 
R is a positive-definite matrix. 
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One may wonder why the extrapolation to infinite 
time is possible in this fashion. As demonstrated be- 
low, linear prediction generates a superposition of oscil- 
lating and exponentially decaying (or growing) terms, 
a type of time dependence that emerges naturally in 
many-body physics: Green's functions of the typical form 
G(k,tu) = (lo — e& — S(fc,w)) _1 are in time-momentum 
representation dominated by the poles; e.g. for a single 
simple pole at lo = u>\ — ir\\ with residue c±, it will read 
G(fc, t) = cie~ luJlt ~ Vlt , and similarly it will be a superpo- 
sition of such terms for more complicated pole structures. 
So the ansatz of the linear prediction is well suited for 
the typical properties of the response quantities we are 
interested in. Note, however, that the method is not gen- 
erally applicable for genuine noncquilibrium situations, 
as e.g. for the evolution of a local observable after a non- 
infinitesimal quench of system parameters. In such cases, 
typically, too many different frequencies can contribute 
to the signal, making the ansatz inappropriate. 

To see the special form of time-series generated by the 
prediction, we introduce vectors x n :~ [x n _i, . . . , x n - p ] T 
such that ([5]) takes the form 



-n+l 



A • x r 



(8) 



with 



A = 



-ai 
1 








-0.2 —a 3 



1 



1 



-dp 









(9) 



Prediction therefore corresponds to applying powers of A 
to the initial vector x^. A (right) eigenvector decompo- 
sition of A with eigenvalues leads to 



XN+r, 



[A 



X N \i = 



(10) 



where coefficients Cj are determined from x^ and the 
eigenvectors of A. The eigenvalues a, encode the physi- 
cal resonance frequencies and dampings. The connection 
is given as on — e ,1JiAt "' iiAt . Excluding exponentially 
growing terms as unphysical, all should obey \ai\ < 1. 
In numerical practice, eigenvalues |aj| > 1 may occur, 
and various remedies exist, all based on the assumption 
that the corresponding a are small. Common manipula- 
tions of occurring > 1 are at — > ai/\ati\, at — > 1/cc*, 
and ai — ► 0. We tested all three and settled on the last 
option. When linear prediction is applicable, the choice 
should not matter (see discussion at the end of Sec. lIII AI . 

From the study of several examples, we found that the 
error of the linear prediction is roughly a function of p-At, 
i.e. p should be adapted to the choice of the time step At 
in the DMRG time evolution. In the simulations, p- At = 
ifit/2 was chosen as a compromise between a large p to 
allow for a superposition of as many different oscillations 



as possible and a small p to have good statistics for the 
fit of the coefficients a,. 

At T = 0, critical one-dimensional systems exhibit 
power-law decays in their time-dependent correlators. 
The superposition of exponential decays is then taken 
to mimic these power-laws 2 ^. At finite temperatures, 
time-dependent correlators S(k, t) decay typically expo- 
nentially for large times (due to thermal broadening), 
making linear prediction especially well-suited for this 
situation. 



III. APPLICATIONS 



XX model 



The Hamiltonian of the XX model reads 



Hxx - ^2(3 fS: 



i+l 



qy qy i 



(ii) 



where we choose the critical field h = — 1. We con- 
sider the observable S(i, t) = ^{[Sj{t), S^(O)}) , which 
is the spectral function of the corresponding model of 
hardcore bosons. This allows for a detailed analysis of 
the method, as correlators for this model can be calcu- 
lated numerically exactly by the evaluation of Pfaffian 
determinant s 26 ! 29 ' 30 . 

The DMRG simulation employed a fourth-order 
Trotter-Suzuki decomposition and a time step At = 



k=it/4, exact 
k=it/2, exact 
k=3it/4, exact 
DMRG 
linear prediction 
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FIG. 2: (Color online.) XX model at /3 = 10: The exact 
Re S(k, t) and the result of DMRG simulation followed by 
linear prediction are in excellent coincidence. The inset shows 
the long-time behavior of the deviation to the exact result. At 
t — tobs, the DMRG simulation was stopped (gray area). For 
the prediction, the deviation (|6]) of predicted to simulated 
values was minimized on the interval (fobs — tfitjtobs] (dark 
gray) with respect to the coefficients ai in ([5]). The hatched 
interval (t a ^ B — p ■ Af, t b a ], determines the value of the first 
predicted value at t = t obs + At. 
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FIG. 3: (Color online.) Lineshapes for the XX model at 
various values of the quasimomentum k (from left to right: 
k = tt/8, 7r/3, and 3tt/4) at f3 = 10 and /3 = 50 : ex- 
act (solid lines) and predicted (dots) lineshapes are in excel- 
lent agreement. For comparison, Parzen-filtered lineshapes 
(dashed lines) obtained from (raw) DMRG data are shown. 



0.2. Deviations | ||p2ff ct ) - \p^f KG )\\/At between ex- 
actly time-evolved state and its DMRG approximation 
were bounded from above by 0.005 per time unit. We 
used a lattice of L = 100 sites 4 ^. To demonstrate the 
potential of the prediction method, the calculations was 
stopped at times i bs = 10.4 for (3 = 10 and i b s = 26.6 
for (3 = 50 (which can be reached with very moderate 
computational resources corresponding to DMRG block 
Hilbcrt spaces of dimension m < 200). For the linear 
prediction we used p = 14 for (3 = 10 and p = 34 
for (3 = 50. The fitting interval was taken from half 
the maximal computation time onward (tat = ^obs/2 and 
p- At = t obs /4). 

In Fig. [5] it can be seen that very accurate predic- 
tions are possible far into the future. As Fig. [3] shows, 
for both temperatures, lineshapes can be calculated very 
accurately upon Fourier transformation of the predicted 
data, allowing for precise experimental analysis. In con- 
trast, even upon windowing (here done by a Parzen fil- 
ter), the raw data provide completely wrong line maxima 
and widths. Almost the whole spectral function is repro- 
duced very exactly. The sole exception occurs for values 
of k very close to zero and 7r: here, the group velocity 
vanishes, dej/dfc — » 0, leading to a very slow decay of 
correlators which is not captured perfectly for the chosen 
temperatures and maximum simulation times t bs- 

In Sec. Ill Bl we gave already arguments for the param- 
eter choice ifit = t t, s /2 and p • At = tgt/2 for the predic- 
tion method. There is one further parameter involved: 
For the determination of prediction coefficients a; in ifS]). 
the autocorrelation matrix R has to be inverted, l[7|). and 
requires some regularization. To this purpose, one may 
cither add a small regularization constant before the in- 
version i?" 1 — > (i?+el) _1 , or project out the cigenspaces 



with eigenvalues below e i.e. i?" 1 — > (P e i?P e ) _1 (pseudo 
inverse) . 

In cases where the linear prediction method is not ap- 
propriate, or the £ Q bs achievable with i-DMRG is too 
small, the results of the prediction are generally sensi- 
tive to variation of the parameters tfn, Pi and e. In this 
respect, the pseudo inverse approach i?" 1 — > (P £ RP e )~ 1 
is favorable, as it produces stronger variations in the re- 
sult for small changes in e, if the linear prediction is not 
applicable. If the method is applicable, the regularization 
parameter e can be varied over several orders of magni- 
tude without effect, but should of course be chosen as 
small as possible (10 _7 -10~ 6 in our calculations). An ef- 
ficient procedure to fix e is to examine plots of the eigen- 
values cti(k) of © as functions of the quasimomentum 
k for several e. For too small e, noisy scatter and many 
abrupt jumps appear, resulting in completely wrong line- 
shapes for some k- values; e is sufficiently big when a^(fc) 
shows smooth "band structures" with only a few abrupt 
jumps. In cases where this is only possible with large e, 
the prediction method is not applicable. 

B. Isotropic Heisenberg spin-i antiferromagnet 

For the Heisenberg antiferromagnet (XXX model) 

^HAFM = Si ■ Si+i (12) 

i 

we have calculated the longitudinal spin structure factor 
(spectral function) S zz (£,t) = (SI(t)S§(O)) . This case 
is interesting and challenging as already at T = a whole 
continuum of excitations contributes to the spin struc- 
ture factor as opposed to the simple dispersion of the XX 
model with a single peak in S(k, u>) for each momentum. 
As concluded from (T = 0) Bethe ansatz calculations, 
the dominant contributions stem from a two-spinon con- 
tinuum bounded from below and above b} r 31 i 32 

= ^| sinfcj and eu(k) = tt\ sin — |. (13) 

The exact contribution of the two-spinon continuum was 
derived in Ref. I331I341 and the contribution of the four- 
spinon continuum computed in Ref.[H (see also Ref. l36h. 

The i-DMRG calculation was done for a lattice of 128 
sites 40 , with a fourth-order Trotter-Suzuki decomposition 
and time steps Aj3 = At = 0.125. The truncation er- 
ror x 2 = J2i> m (X is the 2-norm of the discarded 
Schmidt coefficients Aj) was chosen as \ 2 < 10" 12 during 
the cooling and \ 2 < 10" 10 during the real-time evolu- 
tion, resulting in a maximum DMRG block Hilbert-space 
dimension of to < 1200 for the maximum simulation 
times £ bs = 6.25,7.25,10.5 at (3 = 1,4,16, respectively. 
For the linear prediction ([5]), we used correspondingly 
p = 13, 15, 21 (p • At = Us/4) and t m = t ohs /2. 

The result of DMRG simulation, combined with linear 
prediction is displayed in Figs. 0] and The structure 
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FIG. 4: (Color online.) The longitudinal spin structure factor 
S zz (k, u>) for the isotropic Heisenberg antiferromagnet at j3 = 
4 as obtained from DMRG time evolution followed by linear 
prediction with f b s = 7.25, tfi t = 3.5, and p = 15. Blue lines 
mark the bounds (|13[1 on the (T = 0) two-spinon spectrum. 
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FIG. 5: (Color online.) Line shapes of the longitudinal spin 
structure factor S zz (k,ui) for the isotropic Heisenberg anti- 
ferromagnet as obtained from DMRG and linear prediction 
at finite T, compared to the T = Bethe ansatz result. The 
Bethe ansatz data ("B.A.", provided by J.-S. Cau»2^) com- 
prises the two- and four-spinon contributions, accounting for 
about 98% of the total weight. 



factor converges for (3 — > oo to the T = Bethe ansatz 
results from Ref. 

Fig. [6] illustrates robustness of the linear predic- 
tion against variation of the corresponding parameters. 
Changing the regularization parameter by a factor of 10 
or the length p ■ At of the prediction interval by a factor 
1.5 has no visible effect. 

There is no Bethe ansatz result available for T > 0. 
We hence compare in Figs. [7] and [5] to recent quantum 
Monte Carlo (QMC) data^i. QMC simulations yield the 
spin structure factor S ,lu (k 1 T) on the imaginary-time 
axis. This observable can be determined to very high 
precision. The actual quantity of interest is however 



1 

k=0.8ii 


i 

J \(3=16 


I I I 

P=Po. e " E 
p=1 .5p , c=e 
p=p , e=10e + 














i i i 



0.5 1 1.5 2 

frequency m 



2.5 



FIG. 6: (Color online.) Testing the robustness of the linear 
prediction method at the example of S zz (k,ui). The black 
curves (p = po and e = eo) correspond to the same DMRG 
data as in Fig. [5] (A/3 = At = 0.125; for (3 = 1, A, 16, we used 
po = 13, 15, 21 and e = 10~ 7 , 1CT 7 , 1CT 6 , respectively). The 
blue tilted crosses correspond to simulations where the length 
p- At of the prediction interval was increased by a factor of 1.5. 
The red crosses correspond to calculations where the regular- 
ization parameter e for the inversion of the autocorrelation 
matrix R was increased by a factor of 10. 
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FIG. 7: (Color online.) Comparison of the longitudinal spin 
structure factor S zz {k,uj) for the isotropic Heisenberg antifer- 
romagnet as obtained from DMRG and linear prediction on 
the one hand, and QMC and the maximum entropy method 
on the other hand (provided by S. Grossjohann^ 7 -). The max- 
imum deviation is about 20%. Labels "Bryan" and "his- 
toric" refer to the maximum entropy methods of Refs.l38l and 
Ref. [39l . respectively. 



S^(k,Lu), where 

S^(k,r) = 



1 

2^ 



due-^S^fcu). (14) 



To obtain S^ u {k,uj) from S^ v (k,r), one can employ sev- 
eral variants of the maximum entropy method (see e.g. 
Ref. H for a review). The longitudinal spin structure 
factor S zz (k,oj), Fig. [3 shows a maximum deviation of 
about 20%. If one compares however on the imaginary- 
time axis, S zz (k,T) given in Fig. El QMC and DMRG 
with linear prediction (the DMRG data was transformed 
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FIG. 8: (Color online.) Comparison of the longitudinal spin 
structure factor S zz (k, r) on the imaginary-time axis. DMRG 
with linear prediction (lines) and QMG 3 - 7 - (dots) agree very 
well (maximum discrepancy of ~ 5 ■ 10~ 4 ). It is to be con- 
cluded that the deviations in S zz {k,uj) stem from the ill- 
conditioned analytic continuation which has to be carried out 
for the QMC data and amplifies small discrepancies. 



according to (|14[) ) agree with a maximum deviation of 
~ 5 • 1CU 4 . The discrepancy of the corresponding data for 
S zz (k,oj), is hence to be attributed to the fact that the 
maximum entropy method, employed to transform the 
QMC data to frequency space, is ill-conditioned. This is 
due to the exponential in (|14j) . Small errors of S zz (k,r) 
get blown up by the numerical analytic continuation. 



IV. DISCUSSION 



We have demonstrated how a combination of finite- 
temperature i-DMRG and linear prediction allows for a 
very accurate calculation of spectral functions, at finite 
temperatures. The method is for one-dimensional sys- 
tems favorable to the usual QMC approach: As i-DMRG 
has no sign problem, fermionic systems are also directly 
accessible. Also, no analytic continuation of correlators 
as in QMC calculations is necessary. As opposed to some 
other approaches, the proposed method works over the 
entire temperature regime. An attractive feature of the 
method is provided by the increasing availability of t- 
DMRG codes - which need only minor modification to 
simulate at finite temperatures - and the simplicity of 
the numerically inexpensive prediction method. 
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